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1. Introduction 

Daubechies wavelet matrices Dm are perfect reconstruction orthogonal filter banks 
to which there correspond the orthonormal bases of compactly supported wavelet func- 
tions {2~5^(2 l x— j)}ijez- However, in most of practical applications of these wavelets, 
what matters is the coefficients of Dm and not the form of the corresponding function 
ip. In mere approximation of the irrational coefficients of D m by rational numbers, 
which is desirable in order to simplify the related calculations on a digital computer, 
the perfect reconstruction property of the filter bank Dm (which is its most important 
property) is not preserved in general. Obtained in this way Dm rj Dm is a perfect 
reconstruction orthogonal filter bank only approximately. In the present paper, we de- 
scribe a procedure of approximation of Daubechies wavelet matrices Dm by filter banks 
Dm with rational coefficients which have the perfect reconstruction property exactly. 
This approach depends on a recent parametrization of compact wavelet matrices [B] 
which was developed in parallel with a new matrix spectral factorization method [TJ, 

The paper is organized as follows. The necessary notation and definitions are in- 
troduced in the next section. In Section III, an exact formulation of the problem 
solved is given. In Section IV, the mathematical background of the proposed method 
is provided and the method itself is described in Section V. Some results of numerical 
simulations are presented in Section VI. 

2. Notation and Basic Definitions 

The sets of integer, rational, real and complex numbers are denoted by Z, Q, M, 
and C, respectively. D := {z G C : \z\ < 1} and T := <9B. 5y stands for the Kronecker 
delta, 5ij — 1 if i — j and otherwise, and Idx is the identity map on a set A. 

Let M.p{ m x A) be the set of m x A matrices with entries from a field F (if F is 
omitted it is always assumed that F — C). A row of coefficients (cq, ci, . . . , cm-i) will 
be sometimes called a filter and an m x A matrix will be called an m-channel filter 
bank (with A taps). 

For M e jM(m x A), let M be the matrix with conjugate entries and M* = M . 



1 This method is currently patent pending. 
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U eM(mx m) is called unitary if UU* = U*U — I m where I m stands for the mxm 
identity matrix, and the set of unitary matrices is denoted by U(m). 

V[F] denotes the set of Laurent polynomials with coefficients from a field F, and 
V N [F] := {J2k=-N °kZ k ■ c k E F,k = -N, . . . , N}. If we write just V, the field of 
coefficients will be clear from the context, in most cases V = V[C\. V + C V is the 
set of polynomials (with non-negative powers of z, Ylk=o °k zk P + ) and P C V 
is the set of Laurent polynomials with negative powers of z, Y2k=i c k z ~ k ^ P~ • We 
emphasize that constant functions belong only to V + so that V + fl V~ = 0. Let also 



v% = v ± nv N . 

For power series f(z) = J27=-oo c k zk and N ^ ^ let if( z )]~' [f( z )] + i and 
[f( z )}tn denote, respectively, Efc=-oo c k zk , I2T=o c k zk , Efc=-7V c ^ and Y.k=o c k zk 



and we assume the corresponding functions under these expressions if the convergence 
domains of these power series are known. 

Vim x n) denotes the set of m x n (polynomial) matrices with entries from V, and 
the sets V + (m x n), V^[F](m x n), etc. are defined similarly. The elements of these 
sets P{z) = \pij(z)] are called polynomial matrix functions. When we speak about 
continuous maps between these sets, we mean that they are equipped with a usual 
topology. 

For p(z) = Y,k = _ N c k z k e V, let p(z) = Y,k=-N^k z ~^ a nd for P(z) = [p y (z)] E 
V{m x n) let P{z) = [pij(z)] T G V{n x m). Note that P(z) = (P(z))* when z E T. 

Thus usual relations for adjoint matrices like Pi + Pz{z) = Pi(z) + P2( z ), PiP2( z ) = 
P 2 (z)Pi{z), etc. hold. 

A polynomial matrix function XJ(z) E Vim x m) is called paraunitary if XJ(z)XJ(z) = 
I m for each z E C\{0}, and the set of all paraunitary polynomial matrices is denoted 
by VUim). Note that if U E VU{m), then U(*) E U{m) for each z E T. 

An S E Vj<iim x m) is called positive definite if S(z) is such (XS(z)X* > for 
each ^ X E M(l x m)) for almost every z E T. The polynomial matrix spectral 
factorization theorem (see e.g [3], [5]) asserts that every positive definite S E V?qiyn x 
m) can be factorized as 



where S + E V^im x m) and det S + (z) ^ for each z E D. The representation (PQ) is 
unique in a sense that if S'(z) = S' + (z)S' + (z), then there exists U E U(m) such that 



(1) 



S(z) = S + (z)S + (z), z E C\{0}, 




a 



a 



o 



•mN-l 



'mN-l 



\ 



m— 
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a 



TO — 1 




is said to be a wavelet matrix of ranA; m and genus N, A E WM(m, N, F) (see [U)l 
p. 41]) if the shifted versions of the rows of A by arbitrary multiples of m form an 



orthogonal set, that is 

mN—l 

(3) Yl a k+mrai +ms = mSijSrs, r,s£Z, i,j G {0,1,..., m- 1} 

fc=0 

(it is assumed that a\ = whenever k < or k > mN) and 

mN-l 

(4) 2J a k = m ^io, < i < m. 

k=0 

The conditions ((3]) and 01]) are referred to as the quadratic and linear conditions, 
respectively, defining a wavelet matrix. 

Wavelet matrices with genus 1 (and rank m) are called Haar wavelet matrices, 
H(m,F) := WM(m, 1,F). It can be shown that the first row of any H G H(m,F) 
consists of just Is (see [TOj Lemma 4.4.2]) and if Hi,H 2 G H(m,C), then H\ = 

^ y \ H 2 where U G U (m — 1) (see [TDl Corollary 4.4.3]). Consequently, the only 

real Haar wavelet matrix of rank 2 with determinant 1 is 

(5) ^=(-1 } 
For the matrix ([2]), we consider the Fourier series (the ^-transform) of its rows 

mN-l 

(6) h s {z) = a s k z k , s = 0,1,..., m- 1, 

k=0 

and the corresponding polyphase matrix polynomial A(z) G Vp(m x m) defined by 

JV-l 

(7) A{z) = A k z k 



k=0 



o° ■■■a 

u u m-l 



where A k = [a[ m+j ] G M F (mxm), k = 0, 1, . . . , N—l, e.g. A 



va 



m— 1 _m— 1 



We will heavily use the fact that the quadratic condition (EJ) is equivalent to the con- 
dition on ([7]), 

(8) ^(z)l(z) = ml m , 

i.e. v4(z) is a constant multiplier of paraunitary matrix function. This equivalence can 
be checked by direct computations (see [101 p. 43]). Consequently, A(l) is always a 
Haar wavelet matrix (see [T0~j p. 49]). 

It is well known as well that the quadratic condition ()3]) is also equivalent to the 
following condition on dSJ) (see [TD, p. 96]) 

m—l 

(9) Y h r (zz%)h s {zz%) = m 2 5 rs . 

k=0 
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In signal processing applications, if we split a function (signal) / : Z — > C into m 
parts 

OO j 

(10) / r = -(/' a JC r = 0,1, ...,m- 1, 

777- 

S = — OO 

where aj(-) = a r (- - s), s e Z, (see and (/, a r s ) = J2T=~oo f( k )K( k ), (it is 
assumed that a r (fc) = whenever is outside the range {0, 1, . . . , mN — 1}, so that 
only finitely many products in the above sums differ from 0) which corresponds to the 
filtering by each of the rows a r , r — 0,1, . . . ,m — 1, followed by downsampling with 
rate m, then each of the equivalent conditions Q, (jSj), and (|5J) guarantees that / can 
be reconstructed exactly as follows (see pm Theorem 4.4.23]) 

TO— 1 

(11) /r = £)/r 



/ r ■ 

r=Q 



For this reason, a wavelet matrix is a perfect reconstruction filter bank. 

The linear condition (j3J) implies that a constant signal / : Z — > C emerges from the 
first filter in the representation (11 II) . 

It is said that a wavelet matrix (T5]) has a polynomial-regularity degree d if 

mJV-l 

(12) ^ ^X = 0, p = 0,l,...,d, i = l,2,...,m-l. 

The higher this degree, the more zero coefficients appear in the representation fflOl) of 
a smooth signal /. Note that every wavelet matrix has a polynomial regularity degree 
equal at least to 0. 

3. Formulation of the problem 

The Daubechies wavelet matrix Dn (with 2N taps) is the two-channel filter bank 
with real coefficients 



(13) D 



N 



h \ _ ( a b ai bi ... a N -x b N ^ 
hi J \—bN-i ^AT-i —bN-2 ajv-2 ••• _ bo ao 



which together with quadratic and linear conditions (cf. ([9]) and (HI)) 

(14) \ho(z)\ 2 + \h (-z)\ 2 = 4 when \z\ = 1 
and 

(15) h (l) = 2; h 1 (l) = 

where hj(z) = Yll—o 1 hj[k]z k , j = 0,1, has the polynomial-regularity degree N — 1 

(H2D 

2AT-1 

(16) h[k] ■ k p = for p = 0, l,...,N-l. 

k=0 
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The way of construction of such matrices Dn (computation of coefficients in ( fl3l) ) was 
first established by Daubechies pQ and is described in most books on wavelets [2], [9], 
[TU] . To each matrix there corresponds the Daubechies wavelet ip = ipjy which is a 
supported in [— JV + 1, AT] continuous function of certain smoothness (depending on N) 
such that the system {2~%ip(2 l x — j)}, i,j G Z, forms an orthonormal basis of L 2 (R). 
However the forms of wavelet functions ip^ are mostly of theoretical interest, while 
the numerical values of the coefficients of are very important for applications. 
Since they are irrational numbers in general, during the actual calculations on digital 
computers, these coefficients are quantized and thus D N is approximated by D N . It 
may then happen that satisfies the quadratic condition (fT4"]) only approximately. 
As it has been explained in the preceding section, the quadratic condition on a wavelet 
matrix determines the perfect reconstruction property of a filter bank. 

In the present paper, we propose a method of approximation of D N by D N which 
has rational coefficients and satisfies the quadratic and linear conditions fll4p and ( II 5 p 
exactly. It is obvious that Dn will have the maximal polynomial regularity property 
(I16p only approximately. 



The following theorem, which plays a crucial role in the established method, was 
actually proved in [7]. We present here the simplified proof of this theorem. 

Theorem 1. Let N > 1. For any ^ ip G V^, there exists a unique pair of functions 
a, {3 G such that 



Note that this constant should be positive (for a, /3 ^ 0) since a(z)a(z)+ (3(z)f3(z) = 
\a{z)\ 2 + |/3(^) | 2 for each z G T. 

Proof. It follows from f ll9p that 



4. Mathematical Background of the Method 



(17) 
(18) 



a{z)a(z) + f3(z)f3(z) = 1, z G C\{0} 
«(1) = 1; 0(1) = 0, 



and 



(19) 




Moreover, a and (3 satisfy the condition 

(20) |a(0)| + |/3(0)| >0. 

Lemma 1. Let <\19ty be satisfied for tp G and a, /3 G V^. Then 

(21) a(z)a(z) + (3{z)(3{z) = Const, z G C\{0} 



(22) 




Hence 



aa + (3/3 = $ 2 a - $ x /3 =: $ G V 



■+ 
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and since $ = $ it follows that $ is constant. 



□ 



Proof of Theorem 1. Let <p(z) = Y2k=i"Y kz k ^ e gi yen - We provide a constructive 



'f k Z k 

proof how to find a and 0. First we seek for nontrivial polynomials 

N N 



(23) 



k=0 



M = E 



VkZ 



k=0 



which satisfy ( l22p and hence f[T§]) . If we equate all coefficients of negative powers of z 
of functions <3>i and $2 in ( 122]) to and their Oth coefficients to and 1 respectively, 
then we get the following system of equations in the block matrix form 



(24) 



where 



ex - y = 

QY + X = 1 





/o 


7i 


72 • ■ 


■ 7iV-i 


7aA 




fx Q \ 




(yo\ 








( l \ 




7i 


72 


73 • ■ 


In 







Xi 




Vi 












e = 


72 


73 


74 •• 










x 2 


, Y = 


V2 


, = 





, 1 = 







\7jv 





•■ 





°) 




\xn/ 




\Vn) 




w 




W 



If we substitute the first equation of (j2 
(25) 



Y = 6X 

into the second equation, we get 00 X + X = 1, which is equivalent to 
(26) (00 + W)X = l. 



The system (126]) is nonsingular as is symmetric and 00 = 0*0 is positive definite. 
(Furthermore, all eigenvalues of A := 00 + Jjv+i are grater than or equal to 1, and 
hence ||A _1 || < I as well.) Hence, the coefficients x k and y^, k = 0,1,..., N, in 
fl23|) can be determined from ([26]) and (125]) . and the constructed a and (3 will satisfy 
( |T9~]) . The equation ( |2T]) will be accomplished by Lemma 1 and we can achieve ( IT7j) by 
normalization. 

If now the unitary matrix 



(27) 



U 



is not the identity matrix (note that det U 
and (3 by the equation 

a (3\ fa 
-P a := [-P 



a(l) (3(1) 
-(3(1) 5(1) 

1 by virtue of (H 



we can redefine a 



U~ 



and thus the determined a and (3 will satisfy the conditions (TTTj) - (TT9|) . 

Since the determinant of the product in (TT9T) is 1, we have a(z)§ 2 ( z ) ~ (3(z)&\(z) = 1 
for each z G C (see fl22}). Hence §2U§ holds as well. 



The uniqueness of a pair of polynomials a and ft follows from the uniqueness of 



spectral factorization (see the Introduction) since 
factor of 



1 

<p 1 



a ft 
-ft a 



1 

cp 1 



1 p* 
1 



is the spectral 

□ 



Every process described during the construction of a and ft in the proof of Theorem 
1 is stable under small perturbations of the data, which implies the validity of the 
following 

Corollary 1. Let N > 1, and let JJ : — > x be the map defined according 
to Theorem 1 which assigns a pair of polynomials a and ft to each p 6 . Then \\ 
is a continuous map. 

If we take the coefficients of p rational, then the proof goes through without any 
change and the obtained coefficients of a and ft are rational as well. Thus we have the 
following 

Corollary 2. If <p e V N [Q], then the corresponding polynomials a and ft are from 

Theorem 2. Let N > 1. For any pair of polynomials a,ft£ which satisfy (17) 
and (\20(i there exists a unique <p e such that (\1S\) holds. 

The proof of this theorem is also constructive. 

Proof. Define the function / in a deleted neighborhood of as (see ( BUI) ) 

^ft{z) if a(0)^0 



/(*) 



'0W 



a(z) 



if ft(0)^0 



and let us show that 



satisfies f|T9|) . Indeed, <p G V N (as [f{z)\ = [f(z)] N ) and 

V(z) = f(z) - [f(z)} + 

in a deleted neighborhood of 0. Consider the case a(0) ^ (the case ft(0) ^ can be 
treated analogously). Then 

<pa-ft = (f- [f] + )a -ft = ft- [f] + a -ft = -[f] + a 
and _ 

pft + a = (f- [f) + )ft + a= ^ + aa - [f] + ft = - - [f] + ft, 

a a 

which shows that the functions pa — ft and pft + 5 have removable singularities at 0. 
Since we know that these functions are from V, we conclude that actually they belong 
to V + . Thus ( 1221) and consequently ( |T9l) hold. Observe that 



(2f 



p(z) 





1 








+ ft(z) 




a 


N 



or ip(z) = — 





"1" 


+ 






a(z) 




A 


N 
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so that we need to compute the first N coefficients of l/a(z) or l/ft(z) in order to 
construct <f(z). 

To show the uniqueness of f{z) observe that, by virtue of the Bezout theorem (see 
e.g. [9, Theorem 7.6]), there are a ,(3 G V + such that a(z)a (z) + {3(z)flo(z) = 1. 
Hence, if if G V~ satisfies (122]) . then 

a (z) (ip{z)a{z) - P(z)) + Po(cp(z)P(z) + a(z)) G V + 

and 

(29) if{z) = [a (z)P(z) - f3 (z)a(z)}- 

□ 



As in Theorem 1, the process of construction of ip (see fl28l) ) is stable under small 
perturbations of a and /3. Thus we come to 

Corollary 3. Let N > 1 and Qn C Pjy x "P^ fee the set of pairs (a, (3) which satisfy 
( 77 ) and ffgQl) . T/iere £/ie map ]1 : fijy — )• "P^ defined according to Theorem 2, which 
assigns if G to each (a, (3) G Qn, is continuous. 

We can combine Corollaries 1 and 3 as follows 

Corollary 4. If Q° N c Q n is the set of (a,/3) G f2jv which in addition satisfy (|iflj) . 
£/jen ]j is a continuous one-to-one map from V~ onto Q° N such that Y\ ° 1J = hip- 

Proof. If ip & and (a,/3) G are such that f[T9~j) holds, then ]j(y?) = (a, /3) and 

n(«,/3) = ^ ~' □ 

Corollary 5. For any (a,/3) G Q N , let IK"^) = V 3 IKv 9 ) = Then 
(30) (a,P) = (a/,P')U 

where the matrix U is defined by the equation (W) ) . 
Proof. Because of (JTZJ), f7 G W(2) and = [/*. Thus 

a ^ rr _! ( 5(l)a + 5(l)/3 -/3(l)a + «(1)/3 N 



-P a J \-{-/3{l)a + a{l)/3) a(l)a + 



and (5(1)q + /3(l)/3, /3(l)a + a(l)/3) G Since (see ([H) 

[/-i e p+( 2 x 2), 

we have (a, (3)U~ l = JJ(y>). Thus (a, /J)^ 1 = (a', /?') and flU follows. □ 



1 0\ / a /? 
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5. Description of the Method 



Using the coefficients of ({TBI , we can define 



f 



N-1 



a(z) = —= Y akzk 'i @( z ) 



V2 



k=0 



k=0 



As has been mentioned in Section 2, the quadratic condition (TH| is equivalent to the 
condition for the matrix 



to be paraunitary (see (JTJ) , (JH]) ) , and hence ( 1TT1) holds, while ( fT5i) impies that ^4(1) = 
-^H2, where Hi is the Haar wavelet matrix of rank 2 defined by ([5]). 

Since (a(z), {3(z)) G fijv_i, we can construct <p = Y[(, a iP) according to Theorem 2, 
and 



because of Corollary 5. If we approximate ip by pq G V N _ 1 [Q], <p ~ <pq, and construct 
\J((Pq), then U(^q) e ^at-i(Q) b y Corollary 2, and [J(^) w LK^q) b Y Corollary 1. 
Thus the coefficients of U(<^q) • #2 will be rational and they will approximate and 
bk, k = 0, 1, . . . , N — 1 (see [311) . In this way, we can approximate (TTSl by a matrix 
with rational coefficients which satisfy ( !T4|) and ( !T5|) exactly. 

The proposed method can be generalized for wavelet matrices ([2]) of any rank m 
since the generalization of the main result Theorem 1 used in the method is valid for 
m-dimensional matrices as well (see [H Theorem 1]) and at least the formula (|29|) for 
obtaining ip can be generalized as well (see jl]). However, not for any m, there exists 
a Haar wavelet matrix of rank m with rational coefficients which would provide the 
generalization of formula f )3ip . Consequently, we can approximate any wavelet matrix 
A by A with rational coefficients for which equivalent quadratic conditions (E]), ([S]) and 
([2]) hold exactly, while (jl]) only approximately. As has been explained in Section 2, the 
quadratic condition on a filter bank is decisive for it to have the perfect reconstruction 
property. 



To construct explicitly the fractions which are close to coefficients of Daubechies 
wavelet matrices (see (FIS|) ), a program was written in Mathematica 8. A com- 
plete screening of all possible options has been performed in order to select the fractions 
with minimal denominator in the given range. On a 2GHz Intel Core 2 Duo system 
with 2GB RAM running Ubuntu 11 the calculations took less than a second. As it 
was explained in preceding sections, constructed approximate filter banks Dn have the 
perfect reconstruction property. The results of different approximate computations of 
the coefficients of Daubechies scaling vectors (the first rows of Dn) for genus N = 2 
and N = 3 are presented in the tables below. The pth moments of these coefficients, 




(31) 




6. Computer Simulations and Results 
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M p = Yl'k^o 1 h\[k]k p (see (fl~6l) ). which are not exactly anymore because of approxi- 
mation, are also computed and located in the table. These tables are presented only 
for illustrative purposes and the interested readers can produce the different rational 
approximations which might be more suitable for their specific reasons. 
Table 1. N = 2 





k = 




0.683012701892219 


ifw 0.70588 


if « 0.68597 


i§§2|« 0.68348 




k = 1 




1.18301270189222 


1.17647 


§i« 1-18221 


332288 ^ ii qoqq 
280913 ~ -L-J-O^OO 




k = 2 




0.316987298107781 


^« 0.29411 




88913 ^ rioi eci 
280913 ~ u "3 iD0i 




k = 3 




-0.183012701892219 


-i«-0.17647 


-^T«-0.18221 


280913 ~ 0.18288 




Mi « 




0.0 


0.59 


0.008 


0.001 



Table 2. N = 3 





= 




0.470467207784164 


§§§« 0.5502 


§±§f§« 0.48036 


2677170944 _ q .gn/M 
5703228401 UJ 




fe = 1 




1.14111691583144 


1.1324 


5059904 ^ -i -i QnfiO 
4439725 ~ laJ8DO 


6509075712 ^ 1 1 4 1 on 
5703228401 ~ J-.J-IJ-^J 




fe = 2 




0.650365000526232 


li« 0.5913 


572096 _ 611°9 
887945 u.o^izy 


3712561536 ^ (1 65005 
5703228401 u.uouoj 




fe = 3 




-0.190934415568327 


1056 r, 9Q-I 1 
5249 ~ u - zuii 


170688 q i nooo 
887945 U.J.»A" 


1088205184^ r, i nn°n 
5703228401 ~ U.±»UOU 




fe = 4 




-0.120832208310396 


-H«-°-l 4 15 


553427 q 19465 
4439725 U, "™ J 


686504079 r, lorioy 
5703228401 ~ U.IZUO/ 




fe = 5 




0.0498174997368838 


1^ 0.0687 


233261 ^ n ncoco 
4439725 U.UOZOO 


282357873 q q 495 q 
5703228401 u.uioju 




Mi « 




0.0 


0.256 


0.0357 


-0.0040 




M 2 « 




0.0 


1.622 


0.2169 


-0.0239 
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